---
title: "Stigma and Stereotypes"
author: "Erica A. Metheney"
date: '2020-07-03'
output:
  pdf_document: default
  word_document: default
  html_document:
    df_print: paged
classoption: landscape
always_allow_html: yes
---


```{r setup, include=FALSE}
require("knitr")
require("summarytools")
knitr::opts_chunk$set(echo = FALSE, warning=FALSE, message = FALSE, comment=NA,prompt = FALSE, cache = FALSE, results = 'asis')
opts_knit$set(root.dir = "C:/Users/xmeter/ShareFile/Shared Folders/ZambiaCovidSurvey/Data")

st_options(bootstrap.css     = FALSE,       # Already part of the theme so no need for it
           plain.ascii       = FALSE,       # One of the essential settings
           style             = "rmarkdown", # Idem.
           dfSummary.silent  = TRUE,        # Suppresses messages about temporary files
           footnote          = NA,          # Keeping the results minimalistic
           subtitle.emphasis = FALSE)       # For the vignette theme, this gives better results.
                                            # For other themes, using TRUE might be preferable.

st_css()
```

```{r packages}
library(haven)
require(foreign)
require(nnet)
require(ggplot2)
require(reshape2)
#require(kableExtra)
require(tidyverse)
```

```{r import}
Covid_survey <- read_dta("C:/Users/xmeter/Downloads/Stigma/Covid-survey.dta")
```

# Has Covid Outcome Variable

## Examine the Outcome Variable

Recall: 1 = No, 2 = Not Sure, 3 = Refuse to Answer, 4 = Yes

```{r outcome}
prop.table(table(Covid_survey$HasCovid_raw))
dat = subset(Covid_survey,HasCovid_raw != 3)
prop.table(table(dat$HasCovid_raw))
```

Each category is well populated and the variable is nominal, so we can make use of a multinomial model. 

```{r relevel}
#relevel variables
dat$HasCovid_raw = as.character(dat$HasCovid_raw)
dat$HasCovid_raw[dat$HasCovid_raw=="1"] <- "No"
dat$HasCovid_raw[dat$HasCovid_raw=="2"] <- "Not Sure"
dat$HasCovid_raw[dat$HasCovid_raw=="4"] <- "Yes"
dat$HasCovid_raw = factor(dat$HasCovid_raw)
dat$HasCovid_raw = relevel(dat$HasCovid_raw,ref = "No")

#relevel treatment variables
dat$NeighborAge = as.character(dat$NeighborAge)
dat$NeighborAge[dat$NeighborAge == "1"] <- "25 years old"
dat$NeighborAge[dat$NeighborAge == "2"] <- "60 years old"
dat$NeighborAge = factor(dat$NeighborAge)
dat$NeighborAge = relevel(dat$NeighborAge, ref = "25 years old")

dat$NeighborGender = as.character(dat$NeighborGender)
dat$NeighborGender[dat$NeighborGender == "1"] = "Male"
dat$NeighborGender[dat$NeighborGender == "2"] = "Female"
dat$NeighborGender = factor(dat$NeighborGender)
dat$NeighborGender = relevel(dat$NeighborGender, ref = "Male")

dat$NeighborTimeInCommunity = as.character(dat$NeighborTimeInCommunity)
dat$NeighborTimeInCommunity[dat$NeighborTimeInCommunity == "1"] = "Few Months"
dat$NeighborTimeInCommunity[dat$NeighborTimeInCommunity == "2"] = "Many Years"
dat$NeighborTimeInCommunity = factor(dat$NeighborTimeInCommunity)
dat$NeighborTimeInCommunity = relevel(dat$NeighborTimeInCommunity, ref = "Many Years")

dat$NeighborFamilyOrigin = as.character(dat$NeighborFamilyOrigin)
dat$NeighborFamilyOrigin[dat$NeighborFamilyOrigin == "1"] = "Malawian"
dat$NeighborFamilyOrigin[dat$NeighborFamilyOrigin == "2"] = "Mmwenye"
dat$NeighborFamilyOrigin[dat$NeighborFamilyOrigin == "3"] = "Zambian"
dat$NeighborFamilyOrigin = factor(dat$NeighborFamilyOrigin)
dat$NeighborFamilyOrigin = relevel(dat$NeighborFamilyOrigin, ref = "Malawian")

dat$NeighborSymptoms = as.character(dat$NeighborSymptoms)
dat$NeighborSymptoms[dat$NeighborSymptoms == "1"] = "Cough and High Fever"
dat$NeighborSymptoms[dat$NeighborSymptoms == "2"] = "Injured Leg"
dat$NeighborSymptoms[dat$NeighborSymptoms == "3"] = "High Fever"
dat$NeighborSymptoms = factor(dat$NeighborSymptoms,levels = c("Injured Leg","High Fever","Cough and High Fever"))
dat$NeighborSymptoms = relevel(dat$NeighborSymptoms, ref = "Injured Leg")



#other variables
dat$OutsiderTreatment = as.character(dat$OutsiderTreatment)
dat$OutsiderTreatment[dat$OutsiderTreatment=="0"]<- "No Outsider"
dat$OutsiderTreatment[dat$OutsiderTreatment=="1"]<- "Outsider"
dat$OutsiderTreatment = factor(dat$OutsiderTreatment)
dat$OutsiderTreatment = relevel(dat$OutsiderTreatment,ref = "No Outsider")

dat$Above65 = as.character(dat$Above65)
dat$Above65[dat$Above65 == "0"]<- "No"
dat$Above65[dat$Above65 == "1"]<- "Yes"
dat$Above65 = factor(dat$Above65)
dat$Above65 = relevel(dat$Above65,ref = "No")

dat$VeryWorriedCovid = as.character(dat$VeryWorriedCovid)
dat$VeryWorriedCovid[dat$VeryWorriedCovid == "0"]<- "Not Very Worried"
dat$VeryWorriedCovid[dat$VeryWorriedCovid == "1"]<- "Very Worried"
dat$VeryWorriedCovid = factor(dat$VeryWorriedCovid)
dat$VeryWorriedCovid = relevel(dat$VeryWorriedCovid,ref = "Not Very Worried")

dat$LongTimeResident= as.character(dat$LongTimeResident)
dat$LongTimeResident[dat$LongTimeResident == "0"]<- "Not Long-time Resident"
dat$LongTimeResident[dat$LongTimeResident == "1"]<- "Long-time Resident"
dat$LongTimeResident = factor(dat$LongTimeResident)
dat$LongTimeResident = relevel(dat$LongTimeResident,ref = "Not Long-time Resident")

dat$ChronicIllnesHH= as.character(dat$ChronicIllnesHH)
dat$ChronicIllnesHH[dat$ChronicIllnesHH == "0"]<- "No HH Member with Chronic Illness"
dat$ChronicIllnesHH[dat$ChronicIllnesHH == "1"]<- "HH Member with Chronic Illness"
dat$ChronicIllnesHH = factor(dat$ChronicIllnesHH)
dat$ChronicIllnesHH = relevel(dat$ChronicIllnesHH,ref="No HH Member with Chronic Illness")

dat$OutsiderBias = as.character(dat$OutsiderBias)
dat$OutsiderBias[dat$OutsiderBias == "0"]<- "No Outsider Bias"
dat$OutsiderBias[dat$OutsiderBias == "1"]<- "Outsider Bias"
dat$OutsiderBias = factor(dat$OutsiderBias)
dat$OutsiderBias = relevel(dat$OutsiderBias,ref = "No Outsider Bias")
```

# Multinomial Model 

I run the exact same models (main effects, interactions, baselines) that Cecilia ran using a linear probability model, but in a multinomial frame work. 

Each multinomial model is comprised of 2 regressions since our outcome variable, HasCovid_raw, has three levels. We choose the level "No" to be the baseline category. Thus we are estimating the following models: 

$$\log\left(\frac{Pr(Yes)}{Pr(No)}\right) = \text{intercept} + \text{main effects} + \text{interactions} + \text{error}$$
and 

$$\log\left(\frac{Pr(Not Sure)}{Pr(No)}\right) = \text{intercept} + \text{main effects} + \text{interactions} + \text{error}$$

Recall that there are currently two definitions for outsider, one using all three categories (Zambian, Malawian, Mmwenye) and the other which using Malawian or Not Malawian. 

## Main Effects Models

### Three Levels Outsider

```{r maineffect1, echo = FALSE, message = FALSE}
fit.main1 = multinom(HasCovid_raw ~ NeighborAge + NeighborGender + NeighborTimeInCommunity + NeighborFamilyOrigin + NeighborSymptoms,data = dat)
results1 = summary(fit.main1)
z <- summary(fit.main1)$coefficients/summary(fit.main1)$standard.errors
pvals.main1 <- (1 - pnorm(abs(z), 0, 1)) * 2
result.df = data.frame(results1$coefficients[1,],pvals.main1[1,],results1$coefficients[2,],pvals.main1[2,],stringsAsFactors = FALSE)
colnames(result.df) = c("Not Sure Coef","P-Value","Yes Coef","P-Value")
kable(result.df)
```

### Binary Outsider

```{r maineffect2, echo = FALSE, message = FALSE}
fit.main2 = multinom(HasCovid_raw ~ NeighborAge + NeighborGender + NeighborTimeInCommunity + OutsiderTreatment + NeighborSymptoms,data = dat)
results2 = summary(fit.main2)
z <- summary(fit.main2)$coefficients/summary(fit.main2)$standard.errors
pvals.main2 <- (1 - pnorm(abs(z), 0, 1)) * 2
result.df = data.frame(results2$coefficients[1,],pvals.main2[1,],results2$coefficients[2,],pvals.main2[2,],stringsAsFactors = FALSE)
colnames(result.df) = c("Not Sure Coef","P-Value","Yes Coef","P-Value")
kable(result.df)
```

## Treatment Interactions

### Three Levels Outsider

```{r int1, echo = FALSE}
fit.interactions1 = multinom(HasCovid_raw~NeighborAge + NeighborGender + NeighborTimeInCommunity + NeighborFamilyOrigin + NeighborSymptoms + NeighborSymptoms*NeighborFamilyOrigin,data = dat)
results1 = summary(fit.interactions1)
z <- summary(fit.interactions1)$coefficients/summary(fit.interactions1)$standard.errors
pvals.interactions1 <- (1 - pnorm(abs(z), 0, 1)) * 2
result.df = data.frame(results1$coefficients[1,],pvals.interactions1[1,],results1$coefficients[2,],pvals.interactions1[2,],stringsAsFactors = FALSE)
colnames(result.df) = c("Not Sure Coef","P-Value","Yes Coef","P-Value")
kable(result.df)
```


### Binary Outsider

```{r int2, echo = FALSE}
fit.interactions2 = multinom(HasCovid_raw~NeighborAge + NeighborGender + NeighborTimeInCommunity + OutsiderTreatment + NeighborSymptoms + NeighborSymptoms*OutsiderTreatment,data = dat)
results2 = summary(fit.interactions2)
z <- summary(fit.interactions2)$coefficients/summary(fit.interactions2)$standard.errors
pvals.interactions2 <- (1 - pnorm(abs(z), 0, 1)) * 2
result.df = data.frame(results2$coefficients[1,],pvals.interactions2[1,],results2$coefficients[2,],pvals.interactions2[2,],stringsAsFactors = FALSE)
colnames(result.df) = c("Not Sure Coef","P-Value","Yes Coef","P-Value")
kable(result.df)
```

## All Interactions

### Three Levels Outsider

```{r allint1}
fit.allint1 = multinom(HasCovid_raw~NeighborAge + NeighborGender + NeighborTimeInCommunity + NeighborFamilyOrigin + NeighborSymptoms 
                       + Above65 + VeryWorriedCovid + ChronicIllnesHH+ LongTimeResident + OutsiderBias+
                         NeighborSymptoms*Above65 + NeighborSymptoms*VeryWorriedCovid + NeighborSymptoms*ChronicIllnesHH +
                         NeighborTimeInCommunity*LongTimeResident + 
                         NeighborFamilyOrigin*OutsiderBias,data = dat)
results1 = summary(fit.allint1)
z <- summary(fit.allint1)$coefficients/summary(fit.allint1)$standard.errors
pvals.allint1 <- (1 - pnorm(abs(z), 0, 1)) * 2
result.df = data.frame(results1$coefficients[1,],pvals.allint1[1,],results1$coefficients[2,],pvals.allint1[2,],stringsAsFactors = FALSE)
colnames(result.df) = c("Not Sure Coef","P-Value","Yes Coef","P-Value")
kable(result.df)
```


### Binary Outsider

```{r allint2}
fit.allint2 = multinom(HasCovid_raw~NeighborAge + NeighborGender + NeighborTimeInCommunity + OutsiderTreatment + NeighborSymptoms 
                       + Above65 + VeryWorriedCovid + ChronicIllnesHH+ LongTimeResident + OutsiderBias+
                         NeighborSymptoms*Above65 + NeighborSymptoms*VeryWorriedCovid + NeighborSymptoms*ChronicIllnesHH +
                         NeighborTimeInCommunity*LongTimeResident + 
                         OutsiderTreatment*OutsiderBias,data = dat)
results2 = summary(fit.allint2)
z <- summary(fit.allint2)$coefficients/summary(fit.allint2)$standard.errors
pvals.allint2 <- (1 - pnorm(abs(z), 0, 1)) * 2
result.df = data.frame(results2$coefficients[1,],pvals.allint2[1,],results2$coefficients[2,],pvals.allint2[2,],stringsAsFactors = FALSE)
colnames(result.df) = c("Not Sure Coef","Not Sure P-Value","Yes Coef","Yes P-Value")
kable(result.df)
```

# Relationship Between Has Covid and Symptoms

The first chi-squared test compares the answer to HasCovid with all three symptoms. 

The second chi-squared test compares the answer to HasCovid with the symptoms variable made binary where the options are Not Covid (Injured Leg) or Possibly Covid (fever or cough/fever). 

```{r chisq}
#we want to check if there is a relationship between has covid and symptoms
dat$HasCovid_Binary = rep("",dim(dat)[1])
dat$HasCovid_Binary[dat$HasCovid_raw=="Not Sure"]<- "Unsure"
dat$HasCovid_Binary[dat$HasCovid_Binary != "Unsure"]<- "Answered"

dat$Symptoms_Binary = rep("",dim(dat)[1])
dat$Symptoms_Binary[dat$NeighborSymptoms=="Injured Leg"]<- "Not Covid"
dat$Symptoms_Binary[dat$Symptoms_Binary != "Not Covid"]<- "Possibly Covid"

kable(table(dat$HasCovid_raw,dat$NeighborSymptoms))
chisq.test(dat$HasCovid_raw,dat$NeighborSymptoms)

kable(table(dat$HasCovid_Binary,dat$Symptoms_Binary))
chisq.test(dat$HasCovid_Binary,dat$Symptoms_Binary)

#test if the proportion of unsure is higher for the Not Covid Symptoms
prop.test(x = c(245,661),n = c(1517,3041),alternative = "less")
```

# Difference of Proportions Tests

```{r}
legMa.dat = subset(dat,NeighborSymptoms == "Injured Leg" & NeighborFamilyOrigin == "Malawian")
legZa.dat = subset(dat,NeighborSymptoms == "Injured Leg" & NeighborFamilyOrigin == "Zambian")
legMm.dat = subset(dat,NeighborSymptoms == "Injured Leg" & NeighborFamilyOrigin == "Mmwenye")
feverMa.dat = subset(dat,NeighborSymptoms == "High Fever" & NeighborFamilyOrigin == "Malawian")
feverZa.dat = subset(dat,NeighborSymptoms == "High Fever" & NeighborFamilyOrigin == "Zambian")
feverMm.dat = subset(dat,NeighborSymptoms == "High Fever" & NeighborFamilyOrigin == "Mmwenye")
coughHighFeverMa.dat = subset(dat,NeighborSymptoms == "Cough and High Fever" & NeighborFamilyOrigin == "Malawian")
coughHighFeverZa.dat = subset(dat,NeighborSymptoms == "Cough and High Fever" & NeighborFamilyOrigin == "Zambian")
coughHighFeverMm.dat = subset(dat,NeighborSymptoms == "Cough and High Fever" & NeighborFamilyOrigin == "Mmwenye")

#function
prop.se <- function(p,n){
  y = sqrt(p*(1-p)/n)
  return(y)
}

#----------------------------------------------------Point Estimates and Standard Errors

# ##HELP
# phat.MaLeg = sum(legMa.dat$Help,na.rm = TRUE)/(length(legMa.dat$Help)-sum(is.na(legMa.dat$Help)))
# print(phat.MaLeg)
# phat.ZaLeg = sum(legZa.dat$Help,na.rm = TRUE)/(length(legZa.dat$Help)-sum(is.na(legZa.dat$Help)))
# print(phat.ZaLeg)
# phat.MmLeg = sum(legMm.dat$Help,na.rm = TRUE)/(length(legMm.dat$Help)-sum(is.na(legMm.dat$Help)))
# print(phat.MmLeg)
# 
# se.MaLeg = prop.se(phat.MaLeg,length(legMa.dat$Help)-sum(is.na(legMa.dat$Help)))
# print(se.MaLeg)
# se.ZaLeg = prop.se(phat.ZaLeg,length(legZa.dat$Help)-sum(is.na(legZa.dat$Help)))
# print(se.ZaLeg)
# se.MmLeg = prop.se(phat.MmLeg,length(legMm.dat$Help)-sum(is.na(legMm.dat$Help)))
# print(se.MmLeg)
# 
# phat.MaFever = sum(feverMa.dat$Help,na.rm = TRUE)/(length(feverMa.dat$Help)-sum(is.na(feverMa.dat$Help)))
# print(phat.MaFever)
# phat.ZaFever = sum(feverZa.dat$Help,na.rm = TRUE)/(length(feverZa.dat$Help)-sum(is.na(feverZa.dat$Help)))
# print(phat.ZaFever)
# phat.MmFever = sum(feverMm.dat$Help,na.rm = TRUE)/(length(feverMm.dat$Help)-sum(is.na(feverMm.dat$Help)))
# print(phat.MmFever)
# 
# se.MaFever = prop.se(phat.MaFever,length(feverMa.dat$Help)-sum(is.na(feverMa.dat$Help)))
# print(se.MaFever)
# se.ZaFever = prop.se(phat.ZaFever,length(feverZa.dat$Help)-sum(is.na(feverZa.dat$Help)))
# print(se.ZaFever)
# se.MmFever = prop.se(phat.MmFever,length(feverMm.dat$Help)-sum(is.na(feverMm.dat$Help)))
# print(se.MmFever)
# 
# phat.MaCoughFever = sum(coughHighFeverMa.dat$Help,na.rm = TRUE)/(length(coughHighFeverMa.dat$Help)-sum(is.na(coughHighFeverMa.dat$Help)))
# print(phat.MaCoughFever)
# phat.ZaCoughFever = sum(coughHighFeverZa.dat$Help,na.rm = TRUE)/(length(coughHighFeverZa.dat$Help)-sum(is.na(coughHighFeverZa.dat$Help)))
# print(phat.ZaCoughFever)
# phat.MmCoughFever = sum(coughHighFeverMm.dat$Help,na.rm = TRUE)/(length(coughHighFeverMm.dat$Help)-sum(is.na(coughHighFeverMm.dat$Help)))
# print(phat.MmCoughFever)
# 
# se.MaCoughFever = prop.se(phat.MaCoughFever,length(coughHighFeverMa.dat$Help)-sum(is.na(coughHighFeverMa.dat$Help)))
# print(se.MaCoughFever)
# se.ZaCoughFever = prop.se(phat.ZaCoughFever,length(coughHighFeverZa.dat$Help)-sum(is.na(coughHighFeverZa.dat$Help)))
# print(se.ZaCoughFever)
# se.MmCoughFever = prop.se(phat.MmCoughFever,length(coughHighFeverMm.dat$Help)-sum(is.na(coughHighFeverMm.dat$Help)))
# print(se.MmCoughFever)
# 
# ##MOVE FREELY
# phat.MaLeg = sum(legMa.dat$MoveFreely,na.rm = TRUE)/(length(legMa.dat$MoveFreely)-sum(is.na(legMa.dat$MoveFreely)))
# print(phat.MaLeg)
# phat.ZaLeg = sum(legZa.dat$MoveFreely,na.rm = TRUE)/(length(legZa.dat$MoveFreely)-sum(is.na(legZa.dat$MoveFreely)))
# print(phat.ZaLeg)
# phat.MmLeg = sum(legMm.dat$MoveFreely,na.rm = TRUE)/(length(legMm.dat$MoveFreely)-sum(is.na(legMm.dat$MoveFreely)))
# print(phat.MmLeg)
# 
# se.MaLeg = prop.se(phat.MaLeg,length(legMa.dat$MoveFreely)-sum(is.na(legMa.dat$MoveFreely)))
# print(se.MaLeg)
# se.ZaLeg = prop.se(phat.ZaLeg,length(legZa.dat$MoveFreely)-sum(is.na(legZa.dat$MoveFreely)))
# print(se.ZaLeg)
# se.MmLeg = prop.se(phat.MmLeg,length(legMm.dat$MoveFreely)-sum(is.na(legMm.dat$MoveFreely)))
# print(se.MmLeg)
# 
# phat.MaFever = sum(feverMa.dat$MoveFreely,na.rm = TRUE)/(length(feverMa.dat$MoveFreely)-sum(is.na(feverMa.dat$MoveFreely)))
# print(phat.MaFever)
# phat.ZaFever = sum(feverZa.dat$MoveFreely,na.rm = TRUE)/(length(feverZa.dat$MoveFreely)-sum(is.na(feverZa.dat$MoveFreely)))
# print(phat.ZaFever)
# phat.MmFever = sum(feverMm.dat$MoveFreely,na.rm = TRUE)/(length(feverMm.dat$MoveFreely)-sum(is.na(feverMm.dat$MoveFreely)))
# print(phat.MmFever)
# 
# se.MaFever = prop.se(phat.MaFever,length(feverMa.dat$MoveFreely)-sum(is.na(feverMa.dat$MoveFreely)))
# print(se.MaFever)
# se.ZaFever = prop.se(phat.ZaFever,length(feverZa.dat$MoveFreely)-sum(is.na(feverZa.dat$MoveFreely)))
# print(se.ZaFever)
# se.MmFever = prop.se(phat.MmFever,length(feverMm.dat$MoveFreely)-sum(is.na(feverMm.dat$MoveFreely)))
# print(se.MmFever)
# 
# phat.MaCoughFever = sum(coughHighFeverMa.dat$MoveFreely,na.rm = TRUE)/(length(coughHighFeverMa.dat$MoveFreely)-sum(is.na(coughHighFeverMa.dat$MoveFreely)))
# print(phat.MaCoughFever)
# phat.ZaCoughFever = sum(coughHighFeverZa.dat$MoveFreely,na.rm = TRUE)/(length(coughHighFeverZa.dat$MoveFreely)-sum(is.na(coughHighFeverZa.dat$MoveFreely)))
# print(phat.ZaCoughFever)
# phat.MmCoughFever = sum(coughHighFeverMm.dat$MoveFreely,na.rm = TRUE)/(length(coughHighFeverMm.dat$MoveFreely)-sum(is.na(coughHighFeverMm.dat$MoveFreely)))
# print(phat.MmCoughFever)
# 
# se.MaCoughFever = prop.se(phat.MaCoughFever,length(coughHighFeverMa.dat$MoveFreely)-sum(is.na(coughHighFeverMa.dat$MoveFreely)))
# print(se.MaCoughFever)
# se.ZaCoughFever = prop.se(phat.ZaCoughFever,length(coughHighFeverZa.dat$MoveFreely)-sum(is.na(coughHighFeverZa.dat$MoveFreely)))
# print(se.ZaCoughFever)
# se.MmCoughFever = prop.se(phat.MmCoughFever,length(coughHighFeverMm.dat$MoveFreely)-sum(is.na(coughHighFeverMm.dat$MoveFreely)))
# print(se.MmCoughFever)
# 
# ##HAS COVID
# phat.MaLeg = sum(legMa.dat$HasCovid_raw == "Yes",na.rm = TRUE)/(length(legMa.dat$HasCovid_raw)-sum(is.na(legMa.dat$HasCovid_raw)))
# print(phat.MaLeg)
# phat.ZaLeg = sum(legZa.dat$HasCovid_raw == "Yes",na.rm = TRUE)/(length(legZa.dat$HasCovid_raw)-sum(is.na(legZa.dat$HasCovid_raw)))
# print(phat.ZaLeg)
# phat.MmLeg = sum(legMm.dat$HasCovid_raw == "Yes",na.rm = TRUE)/(length(legMm.dat$HasCovid_raw)-sum(is.na(legMm.dat$HasCovid_raw)))
# print(phat.MmLeg)
# 
# se.MaLeg = prop.se(phat.MaLeg,length(legMa.dat$HasCovid_raw)-sum(is.na(legMa.dat$HasCovid_raw)))
# print(se.MaLeg)
# se.ZaLeg = prop.se(phat.ZaLeg,length(legZa.dat$HasCovid_raw)-sum(is.na(legZa.dat$HasCovid_raw)))
# print(se.ZaLeg)
# se.MmLeg = prop.se(phat.MmLeg,length(legMm.dat$HasCovid_raw)-sum(is.na(legMm.dat$HasCovid_raw)))
# print(se.MmLeg)
# 
# phat.MaFever = sum(feverMa.dat$HasCovid_raw == "Yes",na.rm = TRUE)/(length(feverMa.dat$HasCovid_raw)-sum(is.na(feverMa.dat$HasCovid_raw)))
# print(phat.MaFever)
# phat.ZaFever = sum(feverZa.dat$HasCovid_raw == "Yes",na.rm = TRUE)/(length(feverZa.dat$HasCovid_raw)-sum(is.na(feverZa.dat$HasCovid_raw)))
# print(phat.ZaFever)
# phat.MmFever = sum(feverMm.dat$HasCovid_raw == "Yes",na.rm = TRUE)/(length(feverMm.dat$HasCovid_raw)-sum(is.na(feverMm.dat$HasCovid_raw)))
# print(phat.MmFever)
# 
# se.MaFever = prop.se(phat.MaFever,length(feverMa.dat$HasCovid_raw)-sum(is.na(feverMa.dat$HasCovid_raw)))
# print(se.MaFever)
# se.ZaFever = prop.se(phat.ZaFever,length(feverZa.dat$HasCovid_raw)-sum(is.na(feverZa.dat$HasCovid_raw)))
# print(se.ZaFever)
# se.MmFever = prop.se(phat.MmFever,length(feverMm.dat$HasCovid_raw)-sum(is.na(feverMm.dat$HasCovid_raw)))
# print(se.MmFever)
# 
# phat.MaCoughFever = sum(coughHighFeverMa.dat$HasCovid_raw == "Yes",na.rm = TRUE)/(length(coughHighFeverMa.dat$HasCovid_raw)-sum(is.na(coughHighFeverMa.dat$HasCovid_raw)))
# print(phat.MaCoughFever)
# phat.ZaCoughFever = sum(coughHighFeverZa.dat$HasCovid_raw == "Yes",na.rm = TRUE)/(length(coughHighFeverZa.dat$HasCovid_raw)-sum(is.na(coughHighFeverZa.dat$HasCovid_raw)))
# print(phat.ZaCoughFever)
# phat.MmCoughFever = sum(coughHighFeverMm.dat$HasCovid_raw == "Yes",na.rm = TRUE)/(length(coughHighFeverMm.dat$HasCovid_raw)-sum(is.na(coughHighFeverMm.dat$HasCovid_raw)))
# print(phat.MmCoughFever)
# 
# se.MaCoughFever = prop.se(phat.MaCoughFever,length(coughHighFeverMa.dat$HasCovid_raw)-sum(is.na(coughHighFeverMa.dat$HasCovid_raw)))
# print(se.MaCoughFever)
# se.ZaCoughFever = prop.se(phat.ZaCoughFever,length(coughHighFeverZa.dat$HasCovid_raw)-sum(is.na(coughHighFeverZa.dat$HasCovid_raw)))
# print(se.ZaCoughFever)
# se.MmCoughFever = prop.se(phat.MmCoughFever,length(coughHighFeverMm.dat$HasCovid_raw)-sum(is.na(coughHighFeverMm.dat$HasCovid_raw)))
# print(se.MmCoughFever)
```

Ma = Malawian

Za = Zambian

Mm = Mmwenye

## Outcome Variable: Help

### Symptom: Leg

```{r difftest1}
#------------------------------------------------------------------------Pairwise Tests

#--------------HELP

## LEG
MaLegYes = sum(legMa.dat$Help,na.rm = TRUE)
MaLegTotal = length(legMa.dat$Help)-sum(is.na(legMa.dat$Help))
ZaLegYes = sum(legZa.dat$Help,na.rm = TRUE)
ZaLegTotal = length(legZa.dat$Help)-sum(is.na(legZa.dat$Help))
MmLegYes = sum(legMm.dat$Help,na.rm = TRUE)
MmLegTotal = length(legMm.dat$Help)-sum(is.na(legMm.dat$Help))

### Ma v Za

prop.test(x = c(MaLegYes,ZaLegYes),n = c(MaLegTotal,ZaLegTotal),alternative = "greater")

### Za v Mm

prop.test(x = c(ZaLegYes,MmLegYes),n = c(ZaLegTotal,MmLegTotal),alternative = "less")

### Ma v Mm

prop.test(x = c(MaLegYes,MmLegYes),n = c(MaLegTotal,MmLegTotal),alternative = "greater")
```

### Symptom: Fever

```{r difftest2}

## FEVER
MaFeverYes = sum(feverMa.dat$Help,na.rm = TRUE)
MaFeverTotal = length(feverMa.dat$Help)-sum(is.na(feverMa.dat$Help))
ZaFeverYes = sum(feverZa.dat$Help,na.rm = TRUE)
ZaFeverTotal = length(feverZa.dat$Help)-sum(is.na(feverZa.dat$Help))
MmFeverYes = sum(feverMm.dat$Help,na.rm = TRUE)
MmFeverTotal = length(feverMm.dat$Help)-sum(is.na(feverMm.dat$Help))

### Ma v Za

prop.test(x = c(MaFeverYes,ZaFeverYes),n = c(MaFeverTotal,ZaFeverTotal),alternative = "greater")

### Za v Mm

prop.test(x = c(ZaFeverYes,MmFeverYes),n = c(ZaFeverTotal,MmFeverTotal),alternative = "greater")

### Ma v Mm

prop.test(x = c(MaFeverYes,MmFeverYes),n = c(MaFeverTotal,MmFeverTotal),alternative = "greater")
```

### Symptom: Cough and Fever

```{r difftest3}


## COUGH FEVER
MaCoughHighFeverYes = sum(coughHighFeverMa.dat$Help,na.rm = TRUE)
MaCoughHighFeverTotal = length(coughHighFeverMa.dat$Help)-sum(is.na(coughHighFeverMa.dat$Help))
ZaCoughHighFeverYes = sum(coughHighFeverZa.dat$Help,na.rm = TRUE)
ZaCoughHighFeverTotal = length(coughHighFeverZa.dat$Help)-sum(is.na(coughHighFeverZa.dat$Help))
MmCoughHighFeverYes = sum(coughHighFeverMm.dat$Help,na.rm = TRUE)
MmCoughHighFeverTotal = length(coughHighFeverMm.dat$Help)-sum(is.na(coughHighFeverMm.dat$Help))

### Ma v Za

prop.test(x = c(MaCoughHighFeverYes,ZaCoughHighFeverYes),n = c(MaCoughHighFeverTotal,ZaCoughHighFeverTotal),alternative = "greater")

### Za v Mm

prop.test(x = c(ZaCoughHighFeverYes,MmCoughHighFeverYes),n = c(ZaCoughHighFeverTotal,MmCoughHighFeverTotal),alternative = "less")

### Ma v Mm

prop.test(x = c(MaCoughHighFeverYes,MmCoughHighFeverYes),n = c(MaCoughHighFeverTotal,MmCoughHighFeverTotal),alternative = "greater")


```

## Outcome: Move Freely

### Symptom: Leg


```{r difftest4}
#--------------MOVEFREELY

## LEG
MaLegYes = sum(legMa.dat$MoveFreely,na.rm = TRUE)
MaLegTotal = length(legMa.dat$MoveFreely)-sum(is.na(legMa.dat$MoveFreely))
ZaLegYes = sum(legZa.dat$MoveFreely,na.rm = TRUE)
ZaLegTotal = length(legZa.dat$MoveFreely)-sum(is.na(legZa.dat$MoveFreely))
MmLegYes = sum(legMm.dat$MoveFreely,na.rm = TRUE)
MmLegTotal = length(legMm.dat$MoveFreely)-sum(is.na(legMm.dat$MoveFreely))

### Ma v Za

prop.test(x = c(MaLegYes,ZaLegYes),n = c(MaLegTotal,ZaLegTotal),alternative = "greater")

### Za v Mm

prop.test(x = c(ZaLegYes,MmLegYes),n = c(ZaLegTotal,MmLegTotal),alternative = "less")

### Ma v Mm

prop.test(x = c(MaLegYes,MmLegYes),n = c(MaLegTotal,MmLegTotal),alternative = "greater")

```

## Symptom: Fever

```{r difftest5}

## FEVER
MaFeverYes = sum(feverMa.dat$MoveFreely,na.rm = TRUE)
MaFeverTotal = length(feverMa.dat$MoveFreely)-sum(is.na(feverMa.dat$MoveFreely))
ZaFeverYes = sum(feverZa.dat$MoveFreely,na.rm = TRUE)
ZaFeverTotal = length(feverZa.dat$MoveFreely)-sum(is.na(feverZa.dat$MoveFreely))
MmFeverYes = sum(feverMm.dat$MoveFreely,na.rm = TRUE)
MmFeverTotal = length(feverMm.dat$MoveFreely)-sum(is.na(feverMm.dat$MoveFreely))

### Ma v Za

prop.test(x = c(MaFeverYes,ZaFeverYes),n = c(MaFeverTotal,ZaFeverTotal),alternative = "greater")

### Za v Mm

prop.test(x = c(ZaFeverYes,MmFeverYes),n = c(ZaFeverTotal,MmFeverTotal),alternative = "less")

### Ma v Mm

prop.test(x = c(MaFeverYes,MmFeverYes),n = c(MaFeverTotal,MmFeverTotal),alternative = "greater")

```

### Symptom: Cough and Fever

```{r difftest6}

## COUGH FEVER
MaCoughHighFeverYes = sum(coughHighFeverMa.dat$MoveFreely,na.rm = TRUE)
MaCoughHighFeverTotal = length(coughHighFeverMa.dat$MoveFreely)-sum(is.na(coughHighFeverMa.dat$MoveFreely))
ZaCoughHighFeverYes = sum(coughHighFeverZa.dat$MoveFreely,na.rm = TRUE)
ZaCoughHighFeverTotal = length(coughHighFeverZa.dat$MoveFreely)-sum(is.na(coughHighFeverZa.dat$MoveFreely))
MmCoughHighFeverYes = sum(coughHighFeverMm.dat$MoveFreely,na.rm = TRUE)
MmCoughHighFeverTotal = length(coughHighFeverMm.dat$MoveFreely)-sum(is.na(coughHighFeverMm.dat$MoveFreely))

### Ma v Za

prop.test(x = c(MaCoughHighFeverYes,ZaCoughHighFeverYes),n = c(MaCoughHighFeverTotal,ZaCoughHighFeverTotal),alternative = "greater")

### Za v Mm

prop.test(x = c(ZaCoughHighFeverYes,MmCoughHighFeverYes),n = c(ZaCoughHighFeverTotal,MmCoughHighFeverTotal),alternative = "less")

### Ma v Mm

prop.test(x = c(MaCoughHighFeverYes,MmCoughHighFeverYes),n = c(MaCoughHighFeverTotal,MmCoughHighFeverTotal),alternative = "less")

```

## Outcome: Has Covid

### Symptom: Leg

```{r difftest7}
#--------------HAS COVID

## LEG
MaLegYes = sum(legMa.dat$HasCovid_raw == "Yes",na.rm = TRUE)
MaLegTotal = length(legMa.dat$HasCovid_raw)-sum(is.na(legMa.dat$HasCovid_raw))
ZaLegYes = sum(legZa.dat$HasCovid_raw == "Yes",na.rm = TRUE)
ZaLegTotal = length(legZa.dat$HasCovid_raw)-sum(is.na(legZa.dat$HasCovid_raw))
MmLegYes = sum(legMm.dat$HasCovid_raw == "Yes",na.rm = TRUE)
MmLegTotal = length(legMm.dat$HasCovid_raw)-sum(is.na(legMm.dat$HasCovid_raw))

### Ma v Za

prop.test(x = c(MaLegYes,ZaLegYes),n = c(MaLegTotal,ZaLegTotal),alternative = "less")

### Za v Mm

prop.test(x = c(ZaLegYes,MmLegYes),n = c(ZaLegTotal,MmLegTotal),alternative = "less")

### Ma v Mm

prop.test(x = c(MaLegYes,MmLegYes),n = c(MaLegTotal,MmLegTotal),alternative = "less")

```

### Symptom: Fever

```{r difftest8}
## FEVER
MaFeverYes = sum(feverMa.dat$HasCovid_raw == "Yes",na.rm = TRUE)
MaFeverTotal = length(feverMa.dat$HasCovid_raw)-sum(is.na(feverMa.dat$HasCovid_raw))
ZaFeverYes = sum(feverZa.dat$HasCovid_raw == "Yes",na.rm = TRUE)
ZaFeverTotal = length(feverZa.dat$HasCovid_raw)-sum(is.na(feverZa.dat$HasCovid_raw))
MmFeverYes = sum(feverMm.dat$HasCovid_raw == "Yes",na.rm = TRUE)
MmFeverTotal = length(feverMm.dat$HasCovid_raw)-sum(is.na(feverMm.dat$HasCovid_raw))

### Ma v Za

prop.test(x = c(MaFeverYes,ZaFeverYes),n = c(MaFeverTotal,ZaFeverTotal),alternative = "less")

### Za v Mm

prop.test(x = c(ZaFeverYes,MmFeverYes),n = c(ZaFeverTotal,MmFeverTotal),alternative = "less")

### Ma v Mm

prop.test(x = c(MaFeverYes,MmFeverYes),n = c(MaFeverTotal,MmFeverTotal),alternative = "less")

```

### Symptom: Cough and Fever

```{r difftest9}

## COUGH HIGH FEVER
MaCoughHighFeverYes = sum(coughHighFeverMa.dat$HasCovid_raw == "Yes",na.rm = TRUE)
MaCoughHighFeverTotal = length(coughHighFeverMa.dat$HasCovid_raw)-sum(is.na(coughHighFeverMa.dat$HasCovid_raw))
ZaCoughHighFeverYes = sum(coughHighFeverZa.dat$HasCovid_raw == "Yes",na.rm = TRUE)
ZaCoughHighFeverTotal = length(coughHighFeverZa.dat$HasCovid_raw)-sum(is.na(coughHighFeverZa.dat$HasCovid_raw))
MmCoughHighFeverYes = sum(coughHighFeverMm.dat$HasCovid_raw == "Yes",na.rm = TRUE)
MmCoughHighFeverTotal = length(coughHighFeverMm.dat$HasCovid_raw)-sum(is.na(coughHighFeverMm.dat$HasCovid_raw))

### Ma v Za

prop.test(x = c(MaCoughHighFeverYes,ZaCoughHighFeverYes),n = c(MaCoughHighFeverTotal,ZaCoughHighFeverTotal),alternative = "less")

### Za v Mm

prop.test(x = c(ZaCoughHighFeverYes,MmCoughHighFeverYes),n = c(ZaCoughHighFeverTotal,MmCoughHighFeverTotal),alternative = "less")

### Ma v Mm

prop.test(x = c(MaCoughHighFeverYes,MmCoughHighFeverYes),n = c(MaCoughHighFeverTotal,MmCoughHighFeverTotal),alternative = "less")

```


